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ABSTRACT 

We present the first three-dimensional (3D) simulations of the large-scale mixing that takes place in the 
shock-heated stellar layers ejected in the explosion of a 15. 5Mq blue supergiant star. The blast is initiated 
and powered by neutrino-energy deposition behind the stalled shock by means of choosing sufficiently high 
neutrino luminosities from the contracting, nascent neutron star, whose high-density core is excised and re- 
placed by a retreating inner grid boundary. The outgoing supernova shock is followed beyond its breakout 
from the stellar surface more than two hours after the core collapse. Violent convective overturn in the post- 
shock layer causes the explosion to start with significant large-scale asphericity, which acts as a trigger of the 
growth of Rayleigh-Taylor instabilities at the composition interfaces of the exploding star. Despite the ab- 
sence of a strong Richtmyer-Meshkov instability at the H/He interface, which only a largely deformed shock 
could instigate, deep inward mixing of hydrogen is found as well as fast-moving, metal-rich clumps penetrat- 
ing with high velocities far into the hydrogen envelope of the star as observed, for example, in the case of 
Supernova 1987A. Also individual clumps containing a sizeable fraction of the ejected iron-group elements 
(up to sever al 1(T 3 M ) are obtained in some models. The metal core of the progenitor is partially turned over 
with nickel-dominated fingers overtaking oxygen-rich bullets and both nickel and oxygen moving well ahead 
of the material from the carbon layer. Comparing with corresponding two-dimensional (axially symmetric; 
2D) calculations, we determine the growth of the Rayleigh-Taylor fingers to be faster, the deceleration of the 
dense metal-carrying clumps in the helium and hydrogen layers to be reduced, the asymptotic clump velocities 
in the hydrogen shell to be higher (up to ^4500 km s -1 for the considered progenitor and an explosion energy of 
10 51 ergs, instead of <2000kms -1 in 2D), and the outward radial mixing of heavy elements and inward mixing 
of hydrogen to be more efficient in 3D than in 2D. We present a simple argument that explains these results as 
a consequence of the different action of drag forces on moving objects in the two geometries. 
Subject headings: hydrodynamics — instabilities — shock waves — supernovae: general 



1, INTRODUCTION 

Besides an experimental confirmation of a core col- 
lapse event throug h the detection of supernova neutrinos 
( Hirata et aQ 119871) . the second most important insight pro- 
vided by observations of Supernova 1987 A was the occur- 
rence of large-scale non-radial flow and extensive mixing of 
chemical species in the envelope of the p rogenitor star during 
the explosion (see, e.g. jArnett et al.lll989b . While SN 1987A 
still remains the most prominent and thoroughly-observed ex- 
ample, observations of many other core-collapse supernovae 
(SNe) have meanwhile provided ample evidence that large- 
scale extensive mixing seems to occur generically in such 
events (see, e.g.. IWang & Whe eler 2008). In particular, the 
modeling of the light curve of SN 1987 A using ID radia- 
tion hydrodynamic calculations requires a large amount of 
mixing of Ni outwar d to the H-He interface and of H in- 
ward i nto the He-core (jWooslevlll98 8; Shigevama & Nomotol 
Il990t iBlinnikov et alJEOOOt lUtrobinl 120041) . Moreover, the 
asymmetry of iron and nickel lines in SN 1987 can be ex- 
plain ed, if Ni is concentrated into many h igh-velocity bul- 
lets (iSpvromilio et all Il990t iLi et al.l 119931) . On the other 
hand, the Bochum event, i.e., the sudden development of 
fine-struct ure in the H a line ab out two weeks after the ex- 
plosion (Hanus chik et al.l Il988l) . implies that a high veloc- 
ity (~ 4700km/s) clump of 56 Ni with a mas s of 10~ 3 M(p 
was ejected into the far hemisphere of SN 1987 ( Utrobin et al. 
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fT995l) . 

Observations of near-infrared He I lines from Type II SNe 
between 50 and 100 days after core collapse imply mix- 
ing of 56 Ni into the h ydrogen envelope dFassia et al.l Il998t 
Fassia & Meikle 1999). These authors further conclude that 
those lines are formed in a clumpy environment. Dense 
knots, indications of ejecta shrapnels, and fila ments seen 
in su p ernova remnants b y HST in the visible (Bl air et all 
l2000t iFesen et al.l l2006al) and by ROSAT, Chandra and 
XMM in X-ravs dAschenbach et al.ll995HHughes et alj2000t 
Miceli et al] 120081) also provide strong evidence that mix- 
ing and perhaps even fragmentation is a common process in 
supernova explosions. Spectropolarimetric observations of 
Type II-P and Type Iln SNe at late epochs, when the hydrogen 
envelope starts thinning, reveal strong evidence fo r a globally 
highly aspherical dis t ribution of th e inner ejecta (Fas sia et al.l 
fl998t Leonard et all l200ll |2006b . Similar results are ob- 
tained _from_s_ge£fropojarimetric_o of Typelb/c 
SNe dKawabata et all 12001 IWang et all 12001 iMaund et all 
2007; Mo diaz et alj|2008t IWang & Wheeler! 120081) . Finally, 
studying the neb ula spectra of SN Ic 2 002ap by the means of 
synthetic spectra, Mazzal iet al.l ( 120071) found evidence for an 
oxygen-rich inner core and 56 Ni at high velocities, suggesting 
a highly aspherical explosion especially in the inner parts. 

Two-dimensional hydrodynamic simulations of non-radial 
flow and mixing in the stellar envelopes of core collapse su- 
pernova progenitors have been performed by several groups. 
The first simulations were started by artificially seeding 
Rayleigh-Taylor instabilities (RTIs) in the mantle and en- 
velope of the progenitor and following their evolution until 
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shock breakout from the stellar surface dArnett et all J 19891 : 
Frvxell et al.ll99lHMuller et al.ll991blatlHachisu et aljl99d 
Herant&Benzl I199U lHachisu et al.1 1l992t iHerant & Benzl 
1992). More recent simulations consistently connect the seed 
asymmetries arising from convective flow in the neutrino- 
heated bubble, and by the standing accretion shock instabil- 
ity (SASI) in the "supernova engine" during th e first sec- 
ond of the explosion (see, e.g.. lJanka et al.| [2007) to the later 
Rayleigh-Taylor and Richtmyer-Meshkov instabilities after 
shock passage through the outer stellar laye rs with application 
to SN 1987A (iKifonidis et alJl200l l2006h . Rayleigh-Taylor 
induced mixing and the amount of fallback occurring dur- 
ing artificially triggered (piston model), and initially spher- 
ically symmetric supernova explosions of ze ro- and solar- 
metall icity 15 and 25 M Q stars were studied by Jogge rsTet alj 
(2009). In all models considered by them the velocities of the 
heavy elements we re less than those observed in SN 1987 A. 
Couc h et alj (12009) explored the observational characteristics 
of jet-driven supernovae in a red supergiant progenitor inject- 
ing 10 5I ergs of kinetic energy and varying the characteris- 
tics of the assumed bipolar outflows (jets). Their simulations 
show the development of fluid instabilities that produce he- 
lium clumps in the hydrogen envelope, and nickel-rich jets 
that may account for the non-spherical excitation and sub- 
structure of spectral lines. 

Only a few attempts to 3D si mulations have been made up 
to now. Nagasa wa et al.l (1988) simulated the adiabatic point 
explosion of an n = 3 polytrope using a 3D smoothed parti- 
cle hydrodynamics (SPH) code, and claimed to have found 
pronounced Rayleigh-Taylor instabilities in that event. How- 
ev er, this claim was n either confirmed by the 3D simulations 
of lMuller etaTI(fT989h using a 3D hydrodynamics code based 
on a h igh-resolution shock-captur ing finite-volume method, 
nor bv lBenz & Thielemann! (119901) who also used a 3D SPH 
code. Both studies found only a weak instability, i.e., no ex- 
tensive clumping and mixing, consistent with linear stabil- 
ity anal ysis . Subsequently, 3D simulations by Miill er et al.1 
( 1 99 1 a) and lYamada eTal ( 1 990) using relatively coarse reso- 
lution and/or considering only a wedge-shaped fraction of the 
star showed that seed perturbations grow strongly at the He- 
H interface when a realistic presupe rnova star model i s used 
instead of an n = 3 polytrope. Later iKane et al.l {2000) stud- 
ied the difference in the growth of 2D versus 3D single-mode 
perturbations at the He-H and O-He interfaces of SN 1987 A. 
They found that the 3D single-mode perturbation grows about 
30% faster than a corresponding 2D one. Although the dif- 
ference between 2D and 3D predicted by their simulations is 
not enough to account for the difference between observed 
56 Ni velocities in SN 1987 A and the results of previous 2D 
simulations of SN 1987 A, their results suggest that non-radial 
flow and mixing in supernova envelopes in full 3D is notice- 
ably different from that predicted by 2D simulations. This 
finding is also supported by laser experime nts and compre- 
hensive Rayleigh-Taylor growth sim ulations ( Anuchi na et al.l 
l2004tlCabod2006tlRemington et al.l2006h . The on ly other 3D 
simulations w e are aware of are those published by lNoro et al.1 
(120021 I2004T) in two short conference proceedings. These au- 
thors simulated the propagation of a ID shock wave resulting 
from an artificially triggered spherical explosion of 10 51 ergs 
until it reaches the He-H interface of a SN 1987A progenitor 
model computed by Shigevama & Nomoto ( 1990). Then they 
perturbed the flow inside the shock radius and followed the 
growth of sinusoidal velocity perturbations with a 3D adap- 




FlG. 1 . — Three-dime nsional, neutrin o-driven explosion simulation at about 
0.5 s after core bounce I Scheck 2007). The model is used as initial data for 
our studies. The outermost, bluish, nearly transparent surface is the super- 
nova shock with an average radius of 1900 km. The bright bumpy surface is 
the interface between rising, high-entropy neutrino-heated matter and lower- 
entropy, accreted post-shock gas. An octant is cut out to show the entropy 
distribution in the rising Rayleigh-Taylor bubbles and cooler downflows (the 
entropies per nucleon vary between about 10 and 21 ]cb from blue to yellow 
and red). The neutron star is visible as a dark grey surface near the center. 

tive mesh hydro-code on a grid with an effective resolution of 
up to 4096 3 zones. For an initial perturbation amplitude of 
10%, fingers of heavy elements penetrate into the hydrogen 
envelope with large velocities. 

Here we simulate for the first time, starting with 3D mod- 
els for t he be ginning of neutrino-powered explosions from 
IScheckl d2007l) . the evolution of a supernova blast in 3D from 
the first hundreds of milliseconds to a time three hours later 
when the shock has broken out of a blue progenitor star. We 
find that not only the growth of Rayleigh-Taylor instabilities 
is different in 3D as concluded from other works (see above), 
but also that the propagation of "clumps" in 3D is different 
from the 2D case. In particular the latter effect is of large 
quantitative importance, because we demonstrate here that 
it determines the long-time evolution and in the end the ex- 
tent of large-scale radial mixing in core-collapse supernovae. 
Quantitatively meaningful simulations of observable explo- 
sion asymmetries therefore require modeling in three spatial 
dimensions. 

The paper is organized as follows. In Sect. 2 we describe 
the computational setup and the initial model used in our sim- 
ulations. Our results are discussed in Sect. 3 where we partic- 
ularly address the dynamical evolution and the radial mixing 
of chemical elements. In Sect. 4 we present a simple analyt- 
ical model that explains the differences found for the clump 
propagation between 2D and 3D quite well. Finally, we sum- 
marize our findings in Sect. 5, and draw some conclusions. 

2. SIMULATION SETUP 

For our simulations we used the explicit finite-volume Eu- 
lerian multi-fluid hydrodynamics code, PROMETHEUS, de- 
veloped by Bruce Fryxell and Ewald Miiller (iFrvxell et al.l 
I199U iMuller et all 1 1 99 1 blah . The code integrates the equa- 
tions of multidimensional hydrodynamics to second order 
accuracy in space and time using the dimensional splitting 
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approach of IStrangl (1 19681), the PPM reconstruction scheme 
dColella &"W oodwar dlll984|). and a Riem ann solver for real 
gases according to IColella & Glazl (119841) . Multi-fluid flows 
(involving the advection of, e.g., several nuclear species; see 
below) are treated by means of a set of additional continuity 
equations using the Consistent M ulti-fluid Advection (CMA) 
method of iPlewa & Muileri (fl999h . 

2.1. Computational grid 

We use a computational grid in spherical polar coordinates 
consisting of 1200 x 180 x 360 zones in r, 6, and (^-direction. 
The equidistant angular grid has a resolution of 0.935° in 
(^-direction and 1° in (^-direction, respectively, covering the 
whole sphere except for a (double) cone with a half open- 
ing angle of 5.8° around the symmetry axis of the coordinate 
system (i.e., 0.0325 tt < 6 < 0.9675 tt, and < <p < 2tt). Re- 
flective boundary conditions are imposed in ^-direction, and 
periodic ones in 0-direction. The clipping of the computa- 
tional grid around the symmetry axis avoids a too restrictive 
CFL time step condition. We have not observed any numeri- 
cal artifacts as consequences of this grid constraint. 

The radial mesh is logarithmically spaced between a time- 
dependent inner boundary, initially located at a radius of 
200km, and a fixed outer boundary at 3.9 x 10 12 cm. It has 
a maximum resolution of 2 km at the inner boundary, and a 
resolution of 4 x 10 10 cm at the outer one (corresponding to a 
roughly constant relative radial resolution Ar/r w 1%). We 
allow for free outflow at the outer boundary, and impose a 
reflective boundary condition at the inner edge of the radial 
grid. 

During the simulation we move (approximately every 100th 
timestep) the inner radial boundary to larger and larger radii 
to relax the CFL condition. This cutting of the computational 
grid reduces the number of (the logarithmically spaced) radial 
zones from 1200 at the beginning to about 400 towards the 
end of the simulation, but involves only the innermost few 
percent of the initial envelope mass. We convinced ourselves 
by means of 2D test calculations that this removal of mass has 
no effect on the dynamics and mixing occurring at larger radii 
(the cut radius exceeds at no time a few percent of the radius 
of the Rayleigh-Taylor finger tips). 

2.2. Initial model 

The initial model used for our 3D si mulation s is bas ed on 
3D hydrodynamic explosion models of Scheck (2007|), who 
followed the onset and early development of neutrino-driven 
explosions in 3D with the same code, numerical setup, and 
input p hysics as described in all details for ID and 2D simula- 
tions in lScheck et al.l (120061) . These runs were started from ID 
stellar collapse models a few milliseconds after core bounce. 
The explosion was initiated and powered by the energy in- 
put of neutrinos in the postshock layer, where vigorous con- 
vective activity was present and supported the revival of the 
stalled prompt shock. A defined value of the explosion en- 
ergy was obtained by suitably choosing the magnitude of the 
neutrino luminosities that were imposed at the retreating inner 
grid boundary, which replaced the excised, contracting dense 
core of the forming neutron star. 

The stellar progenitor considered here is a spherically sym- 
metric model of a 15.5M Q blue supergiant with a radius of 
4 x 10 12 cm (IWooslev et al.llT988l Woosley, private communi- 
cation; see als o Kifonidis et al. 2003). Since the models of 
IScheckl d2007l) included only the central part of the star out to 



the middle of the C/O shell at r = 1.7 x 10 9 cm, the envelope of 
the progenitor had to be added for the long-time simulations 
presented here. 

The 15. 5M(n blue supergiant stellar progenitor model was 
evolved by lWoosley et al.l (119881) using a nucle ar reaction net- 
work including a large set of nuclei. However. IS check! (120071) 
did not consider nuclear burning, and he included only a small 
number of nuclear species (p, n, a, and 54 Mn representing the 
heavy element fraction) in their simulations to save computa- 
tional costs. For our initial model w e kept the nucle ar compo- 
sition of the 3D explosion model of lScheckl (120071) inside the 
shock radius, redefining 54 Mn into 56 Ni. The nuclear compo- 
sition at radii greater than the shock radius was taken from the 
1 5.5M^ proge nitor model. 

Scheck (2007) simulated one 3D hydrodynamic explosion 
model with 2° angular resolution, which reached a final time 
of 0.58 s after bounce, and another one with 3° angular resolu- 
tion, which covered an evolution time of 1 s after core bounce. 
Although both models were exploited for our study, we re- 
port here only on simul ations using th e former model that 
according to the data of Scheck (2007) has an (unsaturated) 
explosion energy of 0.6 x 10 51 erg. The deformation of the 
shock front at t ~ 0.5 s after core bounce and the asymme- 
try of the entropy distribution behind the shock in the form 
of rising, inflating hot plumes and infalling, narrow, cool 
downflows are shown in Fig. [T] The structure of the explo- 
sion at this time looks very similar to comparable results by 
iFrver & Warrenl d2002l) . Besides following the shock's prop- 
agation through the stellar envelope of this initial model with 
our 3D code, we also performed a 3D simulation where we ar- 
tificially increased the explosion energy of this initial model 
to 1.0 x 10 51 erg by extending the thermal energy input as ac- 
complishable by ongoing neutrino heating. 

2.3. Equation of state 

We used a tabulated EOS (iTimmes & Swestyl l2000h that 
considers contributions of an arbitrarily degenerate and rel- 
ativistic electron-positron gas, of a photon gas, and of a set of 
ideal Boltzmann gases, consisting of the eight nuclear species 
(n, p, a, 12 C, 16 0, 20 Ne, 24 Mg, and 56 Ni) included in our ini- 
tial model. Note that the pre-shock layers contained only a 
small re gion o f silicon at the time the explosion model of 
Scheck (2007) was adopted as initial data for our long-time 
supernova runs. Since nuclear burning was not taken into ac- 
count in our simulations, but it is possible that this silicon will 
be explosively burned to nickel, we lumped the small amount 
of silicon and the iron-group elements in the postshock layer 
together to what we followed as 56 Ni in our hydrodynamic 
calculations. 

2.4. Simulation runs 

We performed a set of thirteen 2D and three 3D simula- 
tions and compared the results obtained in three dimensions 
to those of the corresponding 2D models. In this paper we 
focus our discussion on one 3D model that has an explosion 
energy of 10 51 ergs. The main findings for this particular 3D 
model are qualitatively similar to those we obtained for less 
energetic 3D explosions, but the resulting maximum values 
for the asymptotic clump velocities vary with the explosion 
energy roughly as v dump oc y/E^. 

The 2D simulations reported here were started from differ- 
ent meridional slices of the 3D initial explosion model. This 
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ensures as much similarity of the initial conditions as pos- 
sible and the same numerical resolution for the 2D and 3D 
runs. But the differences between the chosen slices give rise to 
some variation of the initial asymmetry and explosion energy 
among the various 2D models, which manifests itself both in 
different clump structures and velocities, and the amount of 
mixing in the stellar envelope. 

A set of 13 meridional 2D slices was placed at azimuthal 
angles of 4, 36, 76, 108, 128, 148, 180, 220, 228, 252, 
292, 324, and 360 degrees. Three of these slices (at 128, 
220, and 228 degrees) were chosen at locations where the 
growth of Rayleigh-Taylor fingers was observed to be partic- 
ularly strong in 3D. This guaranteed that such regions were 
not missed by our comparison of 2D and 3D calculations. 
It turned out that the largest 3D Rayleigh-Taylor structures 
at late stages develop in the directions of the biggest initial 
buoyant bubbles, expanding towards the bottom edge and up- 
per right corner (at six o'clock and one o'clock, respectively) 
of Fig. CD 

Both in the 3D and 2D simulations we neglected the influ- 
ence of gravity on the motion of the ejecta. While this has 
no important impact on the dynamics of the expanding ejecta, 
the amount of fallback is underestimated. However, that way 
we could avoid the accumulation of mass near the inner (re- 
flecting) radial boundary which would have implied a consid- 
erably more restrictive CFL condition. 

3. RESULTS 
3.1. Dynamical evolution 

The dynamical evolution of the explosion after its launch 
by neutrino energy deposition and the growth of instabilities 
in the co nsidered 1 5 .5 Mp , prog enitor were described in much 
detail bv lKifonidis et alj d2003l) . It was shown there that the 
asymmetries associated with the convective activity develop- 
ing in the postshock region during the neutrino-heating phase 
act as seeds of secondary Rayleigh-Taylor instabilities at the 
composition interfaces of the exploding star. At about 100 s 
dense Rayleigh-Taylor fingers containing the metals (C, O, Si, 
iron-group elements) have grown out of a compressed shell of 
matter left behind by the shock passing through the Si/O and 
(C+0)/He interfaces. These fingers grow quickly in length 
and while extending into the helium shell, fragment into bal- 
listically moving clumps and filaments that propagate faster 
than the expansion of their environment. 

While for a 2D model with explosion energy around 
1.8xl0 5I erg the Si and Ni containing structures still move 
with nearly 4000km s _I (oxygen has velocities up to even 
5000km s _I ) at 300 s, a strong deceleration occurs when the 
metal clumps enter the relatively dense He-shell that forms 
after the shock passage through the He/H interface. At t > 
10,000 s the metal carrying clumps have dissipated their ex- 
cess kinetic energy and propagate with the same speed as the 
helium material in their surroundings. In the presence of a 
hydrogen envelope and the corresponding deceleration of the 
shock as it propagates through the inner regions of the hydro- 
gen layer (in which case the helium "wal l" below the He/H 
interface builds up), Kifo nidis et alj d2003l) could not observe 
any metal clumps that achieve to penetrate into the hydro- 
gen envelope, in contrast to what was observed in the case of 
SN 1987A. 

The 2D calculations performed in the course of the present 
work confirm these findings (when the velocities are appro- 



sion energies of E exp < 10 51 erg considered here compared to 
roughly twice this value employed by Kifonidis et al. 2003). 
The evolution as well as the final results for the mass distri- 
butions of different chemical elements in velocity space and 
mass space are in good qu antitative agreement with the mod- 
els of lKifonidis et alJ (l2003h . although the latter had much su- 
perior numerical resolution because of the use of a sophisti- 
cated adaptive mesh refinement algorithm. We are therefore 
confident that at least with respect to gross (quasi- ID) fea- 
tures, the results presented here are numerically converged, 
even if many of the structural details, including plumes and 
b ullets, may not be nume rically converged. 

I Kifon idis et alj d2006l) . investigating explosions also with 
£exp ~2x 10 51 erg, proposed a cure of the metal-H mixing 
problem by invoking a sufficiently large dipolar or quadrupo- 
lar asymmetry of the beginning explosion with a globally as- 
pherical shock wave. This was motivated by recent hydrody- 
namical studies that show that large-amplitude, low-^ mode 
oscillations (I defining the order of an expansion in spheri- 
cal harmonics) due to the standing accretion shock instability 
(SASI) can be an important ingredient in the explosion mech- 
anism and in the sequence of processes that leads to observed 
asymmetri es of supernova explosions and neutron star kicks 
(see e.g., iBlondin et alJl2 003: Blondin & Mezzacappai 1 2007 ' 
Scheck et al J 12004 l2006h | Scheckll2007t IScheck et al J 12008 



Burrows et aLT T2006. 20071: iMarek & Jan ka 2009). This on 



priately rescaled with 



to account for the lower explo- 



the one hand led to higher initial maximum velocities of the 
metal-rich clumps and thus their faster propagation through 
the He core on a timescale shorter than the build-up of the 
thick He "wall" that prevented their penetration into the hy- 
drogen shell in the older calculations. On the other hand it 
triggered the growth of Richtmyer-Meshkov instability (RMI) 
at the He/H interface due to the deposition of vorticity by the 
shock. This led to strong mixing of hydrogen inward in mass 
space and down to low velocities, an effect that has turned out 
to be necessary to explain the shape of the li ght curve maxi- 
mum in the case of SN 1987A (IUtrobinll2004l) . 

These conclusions, however, were based on 2D models. 
Therefore the question remained whether 3D might yield dif- 
ferences and whether also in 3D one would have to advo- 
cate a pronounced global low-i 7 mode asphericity of the be- 
ginning explosion as an explanation of the high velocities 
(> 3000km s" 1 ) of nickel clumps in the hydrogen shell and of 
the strong inward mixing of hydrogen observed in SN 1987 A. 
The 3D simulations performed in the present work demon- 
strate that this is not the case and the mentioned observational 
features can be accounted for by 3D models even with explo- 
sion energies around 10 51 ergs for the considered progenitor 
star. 

In order to quantify the global shock deformation of our 
3D initial mod e ls com pared to the 2D cases investigated by 
iKifonidis et alJ (120061) . we did not perform an analysis in 
terms of an expansion in spherical harmonics modes. In- 
stead, sizeable differences become obvious already when a tri- 
axial ellipsoid is con structed as env elope of the bumpy shock 
surface obtained by Scheck (2007) such that the major axis 
of the ellipsoid coincides with the direction of the strongest 
shock expansion. While the 3D model used for the calcula- 
tions of the present paper has an axis ratio of 1 .04: 1 .02: 1.00 at 
t = 0.5 8 s after bounce — another 3D explosion run of lScheckl 
(|2007|) yields 1.16:1 .06: 1 00 at t = 1 s — , t he sho cks of the 2D 
configurations studied by Kifonidis et alJ (120061) typically had 
much larger axis ratios of about 1.5:1.0. 
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FIG. 2. — Surfaces of the radially outermost locations with constant mass fractions of ~ 3% for carbon (green), and oxygen (red), and of ~ 7% for nickel 
(blue). The upper two panels show the directional asymmetries from two different viewing directions at 350 s after core bounce when the metal clumps begin to 
enter the helium layer of the star. The lower two panels display the hydrodynamic instabilities at about 9000 s shortly after the supernova shock has broken out of 
the stellar surface. The side length of the upper panels is about 5 X 10 cm, of the lower panels 7.5 X 10 12 cm. 



The deformed shock causes strong Richtmyer-Meshkov in- 
stability (RMI) to develop at the (C+0)/He and He/H inter- 
fa ces in the 2D models d iscussed here. Compared to the runs 
of iKifonidis et alJ (|2006), however, the smaller asphericity of 
the shock leads to a much slower growth of the instability. In 
our 2D runs the corresponding vortices at the He/H boundary 
at 13000 s have reached a size relative to the average interface 
ra dius similar to what w as observed at 3000 s in the models 
of KifonidisetaLl ((2006). We note in passing that apart from 
resolution-dependent fine structure, the results of the latter pa- 
per could be very well reproduced concerning the RMI vor- 
tex shape, size, and locations by simulations performed for 
the same initial conditions but with the grid setup an d reso- 
lution used in the present work (see Hammer (2009)). The 
growth of the RMI observed in average cases of our 2D cal- 
culations with only weak low-^ mode shock deformation — 



in agr eement with analytic growth rate estimates; lHammerl 
(2009) — is much too slow to produce any significant extent 
of hydrogen-helium mixing so that the final outcome of our 
2D mo dels basically confirms the findings of Kifonidis et al. 
(20030- In the 3D models RMI distortions can be seen at 
the (C+0)/He interface and are likely to contribute to the tur- 
bulent mixing of the metal core with the helium shell of the 
exploding star. At the H/He interface, however, where the 
shock is very close to being spherical, no clear RMI activity 
becomes visible before it is penetrated by fast, metal-carrying 
clumps that have been able to pass through the helium layer 

2 However, there is a considerable spread of the 2D results depending on 
the choice of the meridional plane of the 2D slice and the variation of the 
conditions between the different planes. Slices with somewhat larger initial 
shock deformation show stronger late mixing into the hydrogen shell than the 
more typical "average" 2D slices (see Fig. [8). 
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FIG . 3 . — Volume rendered image of the element distribution at about 9000 s 
(corresponding to the lower two panels in Fig. [2] but displayed for yet an- 
other viewing angle). The blue structures contain a dominant mass fraction 
of nickel, the red fingers mostly oxygen, and green is associated with carbon. 
A mixture of nickel and oxygen appears in pink. Note that the green visible in 
Fig. |2]has been subsumed into the white glow due to the contamination with 
other colors as a consequence of the volume rendering and not of true mixing. 
The side length of the displayed volume is about 7.5 X 10 12 cm. It is remark- 
able that the two long Rayleigh-Taylor fingers from this perspective create the 
impression of ajet and an anti-jet, although no jet-creating mechanism was in 
action at the center of the explosion and the long Rayleigh-Taylor structures 
do actually not point into exactly opposite (i.e., antiparallel) directions (cf. 
Fig.0. 

with still high velocities. 

3.2. Radial element mixing 

Fig. |2] displays the development of these fast-moving 
clumps during our 3D explosion run by showing surfaces 
of constant mass fractions of carbon, oxygen, and nickel 
for two different viewing directions and two different times 
(350 s and ^9000 s after core bounce). Figure [3] provides a 
volume-rendered image of the composition distribution at the 
later time, while Fig. |4]gives composition information on cut 
planes through the mixed stellar core and some of the major 
plumes of different types. Finally, Figs.|6]and|7]present nor- 
malized mass distributions of various nuclear species in the 
radial-velocity and enclosed-mass space for our 3D simula- 
tion compared to an "average" 2D result at several times after 
core bounce, and Fig. [8] provides information for the spread 
of the hydrogen and nickel mass distributions in our set of 2D 
runs at the end of the simulations. 

We stress that what we denote as "nickel" here and in the 
following actually includes the contributions of silicon and 
of the iron-group elements from the postshock region. This 
makes sense because the shock front in the initial model has 
nearly reached the oxygen layer (which implies that only a 
very small fraction of the original silicon is left) and because 
we do not account for nuclear burning in our simulations (for 
which reason no silicon is produced lateron); see Sect. 12.31 In 
the 3D model the total mass of free neutrons is 6 x lO^M©, of 
hydrogen 8.1 6 M©, of helium 5.40 M Q , of carbon 0.12M Q , of 
oxygen O.2OM , of neon 0.048 M Q , of magnesium 0.005 M , 
and of "nickel" O.212M . About 1.35M© of the 15.5M Q 
progenitor are thus absorbed into the neutron star. In the 2D 
runs based on selected meridional slices, the ejecta masses of 



the different nuclear species can differ insignificantly from the 
listed numbers. 

Figure |2] shows that at 350 s oxygen and nickel rich, dense 
fingers, some of which display the typical mushroom shape 
of Rayleigh-Taylor structures, penetrate through the outer 
boundary of the carbon layer of the exploding star. The cor- 
responding 12 C surface is defined by a mass fraction of 3% 
(note that also exterior to this surface, in the surrounding he- 
lium and hydrogen layers, a small but non-zero mass fraction 
of carbon is still present and representative of the heavy ele- 
ments that define the metallicity of the progenitor star). The 
most prominent features at that time have been seeded by the 
biggest high-entropy bubbles visible in Fig. [T] and develop to 
the largest mushrooms visible more than 8000 s later. 

The chemical composition of the various kinds of extended 
plumes in this late stage is very interesting. All of them 
contain high mass fractions of hydrogen and helium, which 
account for about 50 to 70 percent of their mass with spa- 
tial variations in individual filaments. The dominant heavy 
species, however, differs between different structures. The 
largest mushrooms carry predominantly nickel and may con- 
tain a narrow spine with oxygen. In contrast, many of the far 
more numerous, less extended and more narrow fingers are 
composed of a mixture of nickel and oxygen. Also clumps 
with a dominant mass fraction of oxygen can be identified. 
More useful for such an analysis than the surface images of 
Fig. |2] is the visualization by volume-rendering techniques 
in Fig. |2 While the surfaces of constant mass fraction can 
be misleading, because their shape depends strongly on the 
chosen value and no information about the interior composi- 
tion is provided, the ray-tracing image shows that the largest 
Rayleigh-Taylor mushrooms carry mostly nickel, the majority 
of the more extended fingers are composed of a mix of ele- 
ments, and only some clumpy, relatively compact, knotty fea- 
tures (clearly visible along the edge of the pastel-colored and 
whitish glowing core in Fig. [2 contain a dominant fraction of 
oxygen. This is confirmed by the planar cuts through two typ- 
ical Rayleigh-Taylor structures in Fig. [4] The upper left panel 
of this figure, in which a 2D cut through the whole core is 
displayed, demonstrates that on the global scale there are ex- 
tended and separated regions where different species are the 
most abundant chemical constituent. In the core the hydrogen 
and helium admixture is much smaller than in the extended 
fingers and plumes. Except in some smaller "bubbles", where 
H+He dominate, heavies contribute 50-80%, in some clumps 
even more than 90%, of the core composition. 

The time sequences of the normalized mass distributions of 
the nuclear species in the radial-velocity and enclosed-mass 
spac43 of Figs. |6] and |2 visualize the gradual acceleration of 
the outer He and H layers of the exploding star by the out- 
going shock front. Simultaneously, the initially fast metals 
of the core penetrate into the outer composition layers and 
are decelerated. Both the deceleration and the mixing depend 
crucially on the hydrodynamic instabilities that develop at the 
composition interfaces, seeded by the initial explosion asym- 
metries. At 350 s the helium "wall" can be discerned as the 
peak of the distribution between 3500 kms" 1 and 4500 kms" 1 
at an enclosed mass between ~2.5M Q and ~4.5M Q . From a 
comparison of the plots for 2D and 3D it is obvious that nickel 

3 The mass distribution of each nuclear species is normalized such that the 
radial integration in velocity or enclosed-mass space gives unity. Thereby the 
mass AM, of an element ( in a velocity or mass bin is divided by the total 
mass Mi of this element on the grid at the time of evaluation. 
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FIG. 4. — Composition information on 2D cuts through the 3D model structure. Blue corresponds to nickel, red to oxygen, green to carbon. Note that for 
better visibility of the dynamical variation, different color ranges were chosen for the mass fractions in the global-view image (top left) and in the local zooms 
(bottom). Upper left: Mass fractions of Ni, O, and C color coded in a cross sectional plane through the mixed core. The cut plane is indicated in the upper right 
panel and crosses the big mushroom that is enlarged in the lower left panel (this feature grows out of the core at the 12:00 o'clock position), but the basic features 
of the core structure do not depend much on the chosen location of the plane through the grid origin. Upper right: Cutting planes through two representative 
Rayleigh-Taylor structures, the big mushroom growing towards 10 o'clock, and the multi-finger feature extending towards 12:30 o'clock. Both features are 
marked by thick yellow arrows. Lower left: Composition in the cut plane through the big Rayleigh-Taylor mushroom. Lower right: Composition in the cut plane 
through the multi-finger region. All extended plumes and fingers carry a mix of heavy elements with different mass fractions and relative abundances, but also 
large amounts of hydrogen and helium (between ~50 and ~70% of the mass with considerable spatial variations inside individual structures). 



and oxygen experience a much stronger deceleration during 
the subsequent evolution and far less efficient mixing into the 
hydrogen envelope in the 2D case (the plane of the merid- 
ional cut considered in Figs.|6]and|7]is displayed in Fig. HQ. 

4 The set of meridional slices of the 3D initial data that we used for 2D 
simulations exhibits variations at a certain level with respect to gas energies 
and asymmetries of the mass distribution. This leads to corresponding dif- 
ferences in the development of nonradial hydrodynamic instabilities in the 
outer stellar layers. The case picked for Figs.|6]and|2]can be considered as 
"average" concerning metal clump velocities and mixins> but some of the 2D 
slices also yield slightly more extreme results (cf. Fig. [8). In particular, the 
set of slicing directions was chosen such that 3D regions with part icularly 
strong growth of Rayleigh-Taylor fingers were not missed (Sect. l2~4V 



Since neon and magnesium are carbon-burning ashes and well 
mixed with oxygen in the oxygen-layer of the progenitor star, 
the distributions of these three elements closely match each 
other during the mixing evolution and can be hardly distin- 
guished on plots. We therefore have combined them into one 
line in Figs. |6] and [7] 

A closer inspection of Figs. [6] and |7]reveals that the metal 
distributions in both 2D and 3D models peak at the location 
of the metal core (enclosed mass M r < 2.5M Q ) with simi- 
lar velocities (around 2000-3 000 km s _1 at 350s and 1000s, 
around 1500-2500 kms" 1 at 2600s, and around lOOOkms" 1 at 
9000s). In 3D, however, the nickel and oxygen distributions 
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FIG. 5. — Plane of the meridional cut at an azimuthal angle of 148 degrees, 
which was the basis of the 2D simulation whose results are plotted in compar- 
ison to 3D in Figs.|6]and|2] Note that we chose a "typical" (average) direction 
for the cut instead of an extreme one, while the variation of the 2D results in 
different directions is indicated in Fig. [8] 

possess much more massive and wider high-velocity tails up 
to velocities of more than 6000 kms" 1 at 350s. In the long run 
this is causal for the large discrepancies of the radial composi- 
tion mixing. Two effects appear to contribute to the formation 
and continued presence of these high-velocity tails: 

(i) In 3D the Rayleigh-Taylor instabilities at the compo- 
sition interfaces (in particular the (C+0)/He boundary) 
grow more rapidly. This leads to the formation of ex- 
tremely fast (mostly nickel and less oxygen) clumps 
with velocities in excess of 5000 kms" 1 at 350s (up- 
per right panel of Fig. |6]). These dense bullets pene- 
trate deep into the hydrogen layer even before the he- 
lium wall has formed (see the upper two right panels in 
Fig-0- The more rapid growth of the Rayleigh-Taylor 
instabilities also seems to be the reason for the efficient 
mixing of hydrogen into the metal core, which is clearly 
different from the 2D case already at ~ 1000 s. 

(ii) A larger fraction of the nickel and oxygen with a veloc- 
ity of ^4000 kms" 1 at 350s, which is in the ballpark of 
the maximum nickel and oxygen velocities of the 2D 
models, experiences less velocity reduction in the 3D 
simulation. This is obviously linked to a different prop- 
agation and deceleration of the dense metal-containing 
clumps as they move quasi-ballistically through the he- 
lium shell and inner part of the hydrogen envelope un- 
der the influence of drag forces exerted by the surround- 
ing medium. While the fastest bullets with a speed 
of 5000-6000 kms" 1 at 350s still move with 3000- 
4500kms _1 at t > 2600s, the clumps with initially 
^OOOkms" 1 are decelerated to 2000-3000 kms" 1 . In 
contrast, in 2D none of the material of the metal core 
has retained velocities of more than 2000 kms" 1 at the 
time when the shock reaches the stellar surface (f ~ 
8000s). We will return to a detailed discussion of this 
very interesting phenomenon in Sect. |4j where we will 
present simple analytic arguments for an explanation. 



In Fig.[8]the normalized mass distributions of hydrogen and 
nickel in the radial velocity space are plotted at t ~ 9000 s 
shortly after the shock breakout from the stellar surface (when 
the expansion is close to become homologous) for our 3D run 
in comparison to the compilation of results obtained for all 
2D simulations with different chosen meridional slices. While 
the maximum hydrogen velocities are very similar in 2D and 
3D and reach up to about 6000km s" 1 , a clearly bigger hy- 
drogen mass attains the lowest velocities (below 1000 kms" 1 ) 
in the 3D case. Correspondingly, a small fraction of the hy- 
drogen (some hundredths of a solar mass) is mixed deep into 
the metal core of the star to regions with enclosed masses of 
M(r) < 2Mq (see also the right panels of Fig. [7}. These find- 
ings are very similar to the results obtained in 2D by Kifoni- 
dis et al. (2006; cf. figures 6, 8, and 9 there), who considered 
cases with largel y asp herical beginning of the explosion as 
discussed in Sect. 13. 11 1. In 3D similarly strong inward mixing 
of hydrogen does not require a very pronounced global asym- 
metry of the early blast. Richtmyer-Meshkov instability at the 
shock-distorted He/H composition interface plays no crucial 
role in our 3D simulations, because the outgoing supernova 
shock is sig nificantly less deforme d than in the 2D models 
discussed by Kifonid is et alj (120061) . see Sect. 13. II 

Even more impressive than the hydrogen differences is the 
large discrepancy of 2D and 3D with respect to the nickel 
distribution in velocity space and the maximum nickel ve- 
locities. While in 2D even the fastest nickel clumps move 
with just 2200 kms" 1 at t ~ 9000s (Fig. [HJ, a large fraction 
of the nickel retains velocities of more than 3000 kms" 1 and 
some even more than 4000 kms" 1 in 3D (Fig. [8]). From Fig. [6] 
we can conclude that similarly large 2D vs. 3D differences 
are obtained for the oxygen distribution in velocity space, 
whereas the hydrogen and helium distributions of the 2D case 
are close to the results of the 3D run (except for the low- 
velocity tail of the hydrogen distribution mentioned above), 
and the carbon distribution exhibits slight differences only 
between about lOOOkms" 1 and 2500kms"'. Interestingly, a 
significantly larger fraction of the nickel than of the oxygen 
propagates with the highest speeds of clumps. This is com- 
patible with the fact that the metal content of the biggest and 
most extended Rayleigh-Taylor structures in Fig. [2] is domi- 
nated by nickel, while oxygen is concentrated mostly in thin, 
shorter cores of the mushrooms (see Fig.|4|. 

It should also be noted that the spread of the maximum 
nickel velocities obtained in the set of 2D runs in Fig. [8] is 
considerably smaller than the corresponding spread of the 2D 
results in the mass space. The maximum Ni velocities vary 
only between about 1700 kms" 1 for average and 2200 kms" 1 
for extreme 2D cases (in which the meridional slices were 
placed at the locations of the biggest nickel mushrooms in 
3D). In the former models nickel is transported to slightly 
more than 5M , whereas in the latter ones some of the nickel 
gets mixed out to an enclosed mass of even 8-9 M . In all the 
cases, however, the mixing of iron-group elements into the 
hydrogen layer is much less strong in 2D than in 3D. 

Figure [7j showing an average 2D run, confirms these find- 
ings. In the latest stages nickel and oxygen can be observed 
at an enclose d mass of about 5M P) (in rough agreement with 
the results of iKifonidis et al.l (|2003), see their fig. 15), which 
is close to the base of the hydrogen envelop^ In typical 2D 

5 The slightly deeper intrusion of metals into hydrogen in our models can 
be understood as a consequence of the lower explosion energy, which leads 
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FIG. 6. — Normalized mass distributions of hydrogen (black), helium (magenta), carbon (green), oxygen plus neon and magnesium (red), and "nickel" (including 
iron-group elements and silicon; blue) versus radial velocity v r for the 3D simulation (right) and a 2D simulation performed for the meridional slice of the 3D 
model indicated in Fig. \5\ (left). From top to bottom, the distributions are given at about 350s, 1000s, 2600s, and 9000s after core bounce. The binning is done 
in intervals of Av r = lOOkms -1 and the distributions AM,/M, with i being the element index are given per unit length of velocity. Note the large differences 
between the 3D and 2D results of the O and Ni distributions at high velocities and of the hydrogen distribution at low velocities. 
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FIG. 7. — Normalized mass distributions of hydrogen (black), helium (magenta), carbon (green), oxygen plus neon and magnesium (red), and "nickel" (including 
iron-group elements and silicon; blue) versus enclosed mass M(r) for the 3D simulation (right) and a 2D simulation performed for the meridional slice of the 3D 
model indicated in Fig.[5](left). From top to bottom, the distributions are given at about 350s, 1000s, 2600s, and 9000s after core bounce. The binning is done in 
intervals of AM(r) = 0.2Mq and the distributions AMj/Mj with i being the element index are given per unit length of mass. Note that in 3D oxygen and nickel 
clumps are able to penetrate through the hydrogen envelope, while they are stuck in the helium layer in the 2D case. 
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FIG. 8. — Normalized mass distributions of hydrogen (upper panel) versus 
radial velocity v,- and of "nickel" (including iron-group elements and sili- 
con) versus radial velocity v r (middle panel) as well as versus enclosed mass 
M(r) (bottom panel). Results are shown for the 3D simulation (solid line) 
compared to the corresponding 2D results (grey region) at about 9000 s. The 
binning is done in intervals of Avr = 100km s _I and of AM(r) = 0.2Mq, 
respectively, and the distributions AM, /Mi with ;' being the element index 
are given per unit length of velocity or mass, respectively. The grey-shaded 
area indicates the significant variation between the thirteen 2D simulations, 
which were started from different meridional slices of the initial 3D explosion 
model. 



cases the metal-containing clu mps thus get stuck in t he dense 
helium wall as pointed out by Kifonidis et al.l (120031) . In con- 
trast, in the 3D run nickel and oxygen are mixed deep into 
the hydrogen layer (and also a small amount of carbon is car- 
ried from the metal core into the hydrogen shell). The fastest 
nickel and oxygen clumps are able to penetrate nearly to the 



surface of the exploding star (in the mass space). 

The two humps that can be discriminated in the mass dis- 
tribution of these nuclear species well ahead of the bulk of 
the metal core in the mass and velocity space at 9000 s, are 
linked to the few, but big and far-reaching structures, and to 
the greater number of smaller fingers that extend more uni- 
formely to all directions in Fig. [2] The bigger bullets con- 
tain about three to four times more nickel than oxygen. They 
carry in total a mass of ^2x 10~ 2 M© with velocities of 3000- 
4500km s" 1 for the nickel and 3000-4000 km s" 1 for the oxy- 
gen. In contrast, the smaller finger-like structures move with 
1500-3000 kms" 1 and contain a total mass of ~4x 10" 2 M Q . 

It is amazing that in our 3D simulations the fastest of the 
metal clumps have achieved to overtake most of the hydro- 
gen at the time the shock breaks out (at roughly 8000s). 
The nickel and oxygen mixing into the hydrogen envelope is 
clearly more efficient than even in the globally asymmetric ex- 
plosions studied by Kifonidis et al. (2006), where at 10000 s 
the metals were seen to be distributed only up to an enclosed 
mass of about 10M Q (see figure 6 there), although the peak 
metal velocities were roughly the same at ~300s. 

We refrain here from drawing far reaching conclusions on 
the observational consequences of this finding, e.g. concern- 
ing the early visibility of X-ray and 7-ray signatures of such 
strong mixing. While we consider the 2D/3D differences for 
the same explosion model and the same and fixed numeri- 
cal resolution of our first explorative study as enlightening, 
we think that more simulations with different explosion en- 
ergies, different progenitor stars, and in particular also with 
finer grid zoning are desirable before the results can be in- 
terpreted quantitatively with respect to observations of super- 
novae. 

4. A SIMPLE MODEL FOR CLUMP PROPAGATION 

In the following we will present a simple analysis that pro- 
vides support for our hypothesis that the big discrepancies 
found in 3D compared to 2D models are to a large extent the 
consequence of geometry-dependent differences of the propa- 
gation and deceleration of objects ("clumps") under the influ- 
ence of drag forces in the stellar medium of the helium layer. 

4.1. Time-independent toy model 

A first qualitative understanding can be obtained by con- 
sidering objects moving in a stellar medium with a time- 
independent density p(r). We assume that the clump mass 
m = p c V c remains constant and the volume V c and density p c 
of the clumps adjust to the local environmental conditions ac- 
cordingly. With a velocity v relative to the surrounding gas the 
equation of motion of the clumps at a radius r can be written 
as 



dv(r) 

i 

df 



1 dv 2 (r) 
— m 

2 dr 



(1) 



where the drag force for high Reynolds numbers (Re = 
pvV c (A c p) ^> 1 with p being the dynamical viscosity and 
A c the area the object presents to the flow) is given by 



F dag (r) = ^CA c (r)p( r )v 2 (r). 



(2) 



The drag coefficient C is a dimensionless parameter of order 
unity, which depends on the shape and the surface structure of 
the moving object. Equation (Q~|i with Eq. (O has an analytic 
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solution for the velocity as a function of distance: 



with r according to 



v 2 (r) 



6XP ("^/ dr'A c (/)p(/)J =e 



(3) 



when vo is the initial velocity of the object at radial position r$. 
The exponential decay of the initial velocity due to the influ- 
ence of the drag force depends on the density of the medium 
through which the object propagates, on the drag coefficient, 
the mass of the object, and the area it presents to the flow. 

In the following, we will discuss how the term in the ex- 
ponent depends on the dimensionality of the simulation. Be- 
cause of the assumed axial symmetry, in 2D all objects have 
a toroidal form. Let us consider in 2D an ideal torus with ra- 
dius r t , whose center has a distance r from the origin of the 
stellar grid and a distance rsin# from the symmetry axis, and 
compare its outward propagation in radial direction with that 
of a sphere with radius r s in 3D also at radial distance r from 
the center. 

These objects, torus and sphere, are assumed to have a fixed 
shape but their size can vary as they move in radial direction 
from the initial position ro to r. In the beginning, the torus 
radius is r t .o, and the spherical blob has a radius r s ,o- In or- 
der to compare the propagation of structures in 3D and 2D as 
in the hydrodynamical simulations, we suppose that the torus 
and sphere have initially the same cross section and thus ra- 
dius: r s .o = r t o- Since we assume the same density inside the 
objects, pio = /9 s .o, but the torus volume V t .o = 27r 2 r t 2 ro sin 9 in 
general differs from the volume of the sphere, V s o = , the 
masses of torus and sphere are also unequal: m t =/m s . Besides 
these different masses, the two objects also present different 
areas to the flow, A t o = 47rr t o?"osin6> and A s q = 7rr 2 . 

Moving in radial direction, the clumps are assumed to main- 
tain pressure equilibrium with the surrounding medium. As 
a consequence, their density and therefore their volume and 
cross section will change. For the torus one has the scal- 
ing relations V t = V t ,o(r t / r tfi ) 2 (r/ r ) andA t =A t) o(r t /r t ,o)(r/ro), 
while for the sphere they are V s = V , s ,o('sAs,o) 3 and A s = 
A s fl(r s /r s fl) 2 . Using this and the equality of the initial density 
and radii of torus and sphere, the arguments in the exponential 
function of Eq. (01 can be compared for the two geometries to 
yield the ratio: 



Mr) 
Mr) 



Q 8 

C s 37T 



dr'p(r') 



r t r 

n,o r 



dr'p(r') 



r S fl 



(4) 

In this equation the radii of the torus, r t , and of the sphere, 
r s , in general depend on the distance r from the grid center, 
if the clump propagates through a stellar medium with a pres- 
sure gradient. It should be noted that the ratio M r )/M r ) is 
nearly independent of both the mass difference and the initial 
difference of the torus and sphere areas perpendicular to the 
radial direction (the numerical factor 8/(37r) is very close to 
unity). Instead, the deceleration of the motion by drag forces 
depends crucially on the evolution of the areas as reflected by 
the radius-dependent functions in the integrands. 

We assume that the pressure in the clump (torus or sphere) 
scales with the density in the clump as P c oc p" ! oc V~ m , and 
that the pressure in the stellar environment follows a radial 
profile P(r) oc f~ n . Using this and the requirement of pressure 
equilibrium, we find that the radii of torus and sphere change 
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Plugging Eq. (0 into Eq. (0), we obtain 



(5) 
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(6) 



The torus area perpendicular to the radial direction grows 
faster than that of the sphere if n < 3m, in which case the 
integral in the numerator increases more rapidly with r than 
the one in the denominator. This difference is a consequence 
of the fact that any radial motion leads to a larger distance of 
the torus from the polar axis of the grid. This stretches the 
torus and causes an adjustment of its cross section according 
to pressure equilibrium. Since in addition to the area effect 
also the drag coefficient of the toroidal object is likely to be 
larger than that of the sphere, C t > C s , the motion of the torus 
is clearly more strongly damped by the drag than that of the 
sphere if n < 3m. 

In the supernova the clumps propagate through the medium 
of the dense shell of shock-accelerated progenitor matter. In 
this shell the density and pressu re gradients are nearly flat, i.e., 
p ~ const and n ~ (see, e.g.. lNomoto et al] ( 1 19941) . figs. 31 
and 32). This case is particularly suitable to illustrate the in- 
adequacy of axisymmetric toroidal clumps in 2D instead of 
spherical clumps in 3D. While the latter propagate through 
the homogeneous medium without any change of their area 
perpendicular to the direction of the motion, the tori are un- 
avoidably stretched as they reach larger radii r, which leads 
to an increase of the cross-sectional area they present to the 
flow. Equation (0 yields: 
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(7) 



which is larger than unity for r > ro. A stronger deceleration 
of the toroidal clump is thus caused by the artificial constraints 
associated with the axisymmetry of all structures in the 2D 
geometry. 

4.2. Time -dependent toy model 

In a second step we intend to obtain a more quantitative 
insight by considering the clump propagation in the dense he- 
lium layer that forms after the passage of the outgoing shock 
through the He/H interface (for details, see Kifonidis et al. 
2003). The interaction with this shell was identified as the 
main reason why metal-carrying rising Rayleigh-Taylor struc- 
tures were unable to penet r ate int o the H-shell in the 2D sim- 
ulations of Kifonidis et al. (2003). Here we will integrate the 
equation of motion for spherical clumps decelerated by the 
drag force in the 3D environment (in contrast to the 3D case 
of a spherical clump, the equation of motion of a 2D toroidal 
object defies an analytic solution due to the explicit depen- 
dence of the clump volume and cross-sectional area on the 
radial coordinate position r). 

We assume that the He layer is homogeneous but has a 
time-dependent density p(t). This is an acceptable simpli- 
fication because numerical simulations show that the radial 
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density variations in this layer are much smaller than the time- 
dependent decrease of the density due to the expansion of the 
star during the first three hours of the explosion. In the 3D 
case this leaves us just with a time-dependence of the drag 
force, F dra „ = FdragO), and thus 



dv(r) 1 , 
m-±-l =--QA s (t)p(t)v 2 (t). 



(8) 



Again we consider m = const and C s = const and v(f ) as the 
clump velocity relative to the surrounding medium. We fur- 
ther assume that the density of the helium layer evolves like 
pit) = po(k)/t) k and the pressure and density are linked by 
P(t) = Kp" 1 ^). Using the same pressure-density relation for 
the medium inside the clumps, P s = K s p7, but with a different 
constant K s =/K (corresponding to an entropy of the clumps 
different from that of the surroundings), the request of pres- 
sure equilibrium between clumps and environment allows us 
to calculate the time evolution of the clump radius, 



r s (t) = r sfi 



(hp 
Po 



(9) 



and therefore the time-dependent cross-sectional area that 
the clumps present to the surrounding medium: A s (f) = 
A $i0 (r/r S fl) 2 = A*(t /t () ) 2k ^ with constant A*. Inserting this 
and pit) into Eq (8), we can integrate the resulting differen- 
tial equation by separation of variables. This yields the time 
evolution of the velocity (for k ^ 3) as 



v(f) = v Q {\ + 



C S A* 
2m 



■votopo 
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l-k/3 



l-k/3 



1 



(10) 

We specialize now to k = 2 and 7 = f , which is in reason - 
able agreement with the simulations (fKifonidis et al.l [2003). 
Rewriting the constant factor of the second term in the de- 
nominator by using the expressions for clump area and vol- 
ume, we end up with 



v(f) = v <l + 



47rV 
37 



C s vo?o 



v m / 



to 



(11) 

The fastest clumps, which are also the most massive ones with 
m~4x 10~ 3 M Q in our 3D simulations, enter the dense he- 
lium shell with a velocity of vo ~ 2000—4000 km s" 1 (relative 
to the ambient medium) at fo ~ 300s. At this time the density 
of the medium in the shell is po ~ 0. 1 g cm -3 . The clumps typi- 
cally have roughly three times lower entropy than the medium 
in the helium shell so that K s /K ~ 1/3 (most of the numeri- 
cal valu es used here can be ve rified by closely inspecting the 
plots in Kifon idis et al.l (|2003)). Adopting for the drag coef- 
ficient plausible number^ of C s ~ 0.05 . . . 0.2, we obtain at 
t = 9000s velocities of v - ^OO-SlOOkms" 1 for the clump 
motion through the helium material, corresponding to a de- 
celeration by roughly a factor of 2. This means that the fastest 
objects are expected to still speed with ^2400^400 km s -1 
relative to an observer in the lab frame. This is in good agree- 
ment with our hydrodynamical simulations and can thus be 
interpreted as support of our picture of the clump propaga- 

6 Unfortunately, no information is available in the literature for numbers 
that apply to the problem considered here. 



tion under the influence of the drag forc^ In contrast, in 2D 
models the maximum clump velocities are reduced to slightly 
more than lSOOkms" 1 in the lab frame (see Figs. [6j [8] and 
iKifonidis et al.l d2003h ). which is equal to the expansion ve- 
locity of the He-shell matter so that the clumps get stuck in 
their surroundings and are unable to penetrate into the hydro- 
gen envelope (Fig. [7). 

5. SUMMARY AND CONCLUSIONS 

We have presented the first 3D simulations that followed 
the development of mixing instabilities in a supernova over a 
timescale of hours after the explosion was initiated by neu- 
trino energy deposition. The considered progenitor was a 
15.5M Q blue s upergiant star, which had als o been adopted 
for the works of Kifonidi s et al.1 (I2003L 12006). Our long-time 
runs were started a t roughly 0.5 s after core bounce from 3D 
models provided by Scheck ( 2007), who triggered the blast by 
suitably chosen neutrino luminosities imposed at the retreat- 
ing inner grid boundary that replaced the high-density core 
of the nascent, contracting neutron star. The supernova runs 
we focussed our discussion on had an explosion energy of 
10 51 ergs. Our particular goal was a close comparison of the 
growth of mixing instabilities at the composition interfaces of 
the exploding star in 3D and 2D calculations. For the latter, 
we investigated a set of selected meridional slices of the initial 
3D explosion models. 

Our comparison revealed that the asymptotic velocities of 
metal-rich clumps are much higher in 3D than in the corre- 
sponding 2D cases. In the latter the iron-group and oxygen 
carrying dense fingers get stuck in the massive, dense helium 
shell ("wall") that forms below the base of the hydrogen en- 
velope, and stay comoving with the helium matter there at 
a velocity of ^ - lSOO kms" 1 , a fact that was pointed out by 
Kifoni dis et al.1 (120031) . whose results are in agreement with 
our 2D simulations. In contrast, in the 3D runs we obtained 
"nickel" (= elements of the iron group plus siliconf] and oxy- 
gen (plus neon and magnesium) bullets propagating at max- 
imum velocities of 4500kms _I and a large fraction (~10- 
20%) of the metal core of the progenitor star was seen to 
reach velocities of more than 2000 kms" 1 . The most extended 
and fastest Rayleigh-Taylor structures, some of which have a 
mass of several 10~ 3 M Q , were found to contain mostly nickel. 
These nickel "bullets" expand significantly more rapidly than 
the longest fingers with dominant or appreciable oxygen con- 
tent, which are smaller on average but much more numerous 
than the large nickel features. Iron-group nuclei, neon, and 
magnesium are thus carried far into the hydrogen layer. They 
move well ahead of oxygen-rich knots and both iron and oxy- 
gen overtake the material from the carbon layer in the ejecta. 
The onion-shell stratification of the progenitor is thus partially 
turned over during the explosion. 

Besides strong mixing of heavy-elements into the hydrogen 
envelope, we also found about O.O86M of hydrogen being 
transported deep into the metal core of the star to regions with 
velocities of <1000kms _1 , ~0.034M© to an enclosed mass 
of M(r) < 2M Q , and ~O.19M to M(r) < 4M Q . Although 

7 Note that the numerical result for the final velocity varies relatively little 
when different combinations of values for vo, po, and tg within reasonable 
ranges are adopted. 

8 We note again that in our present models without nuclear burning a small 
region of silicon in the collapsing stellar core was lumped together with the 
iron-group nuclei to what we called "nickel". Silicon was therefore not traced 
separately as a nucleosynthetic constituent of the ejecta. 
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this result seems to be roughly compatible with observational 
constraints based on the lightcurve shape of SN 1987 A (see 
Blinnikov 1999 and Kifonidis et al. 2006), we refrain from a 
more detailed comparison and too far reaching conclusions, 
because neither the considered progenitor star is very suitable 
to compare with SN 1987 A, nor were the explosion parame- 
ters (mass cut, nickel mass, explosion energy) adjusted to be 
optimally consistent with the observations of this supernova. 
Moreover, the exact amount of inward hydrogen mixing must 
be expected to be sensitive to the numerical resolution and 
thus at least convergence tests with better spatial zoning are 
re quired. 

Kifo nidis et alj d2006l) obtained similar results on the ba- 
sis of 2D explosion models where a supernova shock with a 
large global asphericity triggered strong Richtmyer-Meshkov 
instability at its passage through the composition interfaces, in 
particular also at the helium/hydrogen boundary. Such an ef- 
fect, however, does not play any significant role in the present 
3D simulations, because the shock starts out with much less 
initial deformation and is close to spherical when it crosses 
from the He into the H shell. Nevertheless, the peak veloci- 
ties of the meta ls in our 3D simulation s are even higher than 
those found by Kifonidis et al.l ( 120061) . and correspondingly 
the fastest nickel and oxygen clumps are able to penetrate 
even farther into the hydrogen envelope. 

A closer comparison showed that a small fraction of the 
metal core receives velocities of more than 6000km s" 1 ini- 
tially, which is significantly higher than in o ur 2D runs (but 
simila r to the peak metal velocities seen by Kifonid is et alJ 
(2006)). These clumps penetrate into the hydrogen envelope 
faster than the helium wall builds up. They achieve to retain a 
speed of more than 4000 km s" 1 at the time of shock breakout 
and e xperience much less d eceleration than in the 2D mod- 
els of [Kifonidis et al] (12006b . Moreover, a larger fraction of 
the core metals start out with velocities of about 4000 km s" 1 , 
which is of the order of the speed of the fastest clumps in our 
2D models, but again the deceleration of these dense bullets is 
less strong than in the 2D case. This suggests that the final dif- 
ferences can be attributed only partly to 2D/3D differences in 
the growth rate of perturbatio ns at the compo sition interfaces, 
which were pointed out by iKane et ail d2000l) . Instead, the 
discrepancies between 3D and 2D results that develop when 
the majority of the metal-rich clumps begins to enter the dense 
helium wall below the He/H boundary after ~ 1000 s, and the 
increase of these discrepancies over the first 1-2 hours after 
the onset of the explosion while the clunps make their way 
through the helium layer, require a different explanation. 

We therefore hypothesized that these differences are a con- 
sequence of differences in the drag forces of the surround- 
ing medium, which affect the propagation of the metal-rich 
clumps and lead to their deceleration as they plow through 
the helium layer. By means of a simple analytic model we 
could demonstrate that the different geometry of the bullets 
in 2D and 3D — toroidal (axisymmetric) structures versus 
quasi-spherical bubbles — can explain the differences ob- 
served in our simulations. Making reasonable assumptions in 
a toy model description of the conditions in the exploding star 
and using plausible values of the drag coefficient, we could 
come up with numbers that are in the ballpark of the asymp- 
totic nickel clump velocities in the hydrogen shell found in 
our 3D hydrodynamic models. 

In this context experimentally and numerically established 
data for the drag coefficient in dependence of the clump ge- 



ometry would have been highly desirable but could not be 
found in the literature. A determination of corresponding val- 
ues may be a valuable goal for future investigations, e.g. at the 
National Ignition Facility (NIF) at Livermore. It must be sus- 
pected that the drag force acting on the dense metal clumps is 
not independent of the grid resolution and insufficient zoning 
could lead to an underestimation of drag effects on the clump 
propagation. Without applying any adaptive mesh refinement 
(AMR) techniques, the 1200 radial zones and 1° resolution 
in azimuthal and 0.935° in lateral direction were the best we 
could allow for as a compromise between accuracy and CPU- 
time requirements. Although our 2D results are in reasonably 
good agreemen t with the 2D AMR calculations performed by 
Kifonidis Hal] (12001 [2006), we still think that simulations 
with varied mesh size or an AMR code would be useful to 
confirm our findings in 3D. 

In this pilot study we were aiming at a first analysis of fun- 
damental differences between mixing instabilities during the 
long-time evolution of supernova explosions in two and three 
dimensions rather than attempting to match the observational 
appearance of special events like, e.g., SN 1987A (although 
we used this supernova as the best studied one for our motiva- 
tion and also referred to it when our findings agreed with ba- 
sic features that SN 1987 A is likely to have in common with 
a wider class of cases). On the basis of our results, which 
appear suitable to account for the strong inward mixing of 
hydrogen and outward mixing of metals that is needed to ex- 
plain the smoothness of the light curve as well as the early 
appearance of X-ray and gamma-ray emission of SN 1987 A, 
we therefore do not conclude that all observational properties 
of this supernova can be reproduced by a relatively weakly 
aspherical initiation of the explosion. More simulations are 
needed to clarify this point, in particular by employing a 
stellar model that is compatible with all knowledge acquired 
about the progenitor structure of SN 1987 A during the past 20 
years. 

While we think that the 2D/3D differences found in our in- 
vestigation are likely to be generic, the extent of the radial 
mixing, the size and mass distribution of the metal fingers, 
their velocities and spatial distribution, and the composition 
of individual clumps must be expected to depend strongly 
on the structure of the progenitor star, the explosion energy, 
and the initial asymmetry of the blast. Having advanced the 
modeling of mixing instabilities in supernova explosions to 
the third dimension, we therefore plan to direct our focus to 
the exploration of a wider variety of progenitor stars with 
the goal to establish a link between explosion models and 
observations. Besides SN 1987A, the geometry of the Cas- 
siopeia A remna nt has been studied observationa l ly in ver y 
much detail (e.g jHwang et all ( 120041) : lFesen et alj d2006albh . 
and references therein). It exhibits properties such as the ra- 
dial turnover of the metal core that are compatible with the 
findings of this work, and a directional asymmetry that can- 
not be accounted for by models constrained to axial symme- 
try. Therefore not only the differences of 3D versus 2D results 
reported in this paper but also observational properties of su- 
pernovae and supernova remnants point to the indispensability 
of modeling the explosions in three dimensions. An interest- 
ing question certainly is, whether the jet-like, silicon-rich fea- 
tures observed in Cas A could be morphologically linked to 
the structures seen as the fastest and most massive Rayleigh- 
Taylor fingers in our simulations. Although our results appear 
suggestive in this respect, a convincing answer requires inves- 
tigations of progenitors and explosion parameters fitting this 



Three-dimensional supernova explosion models 



15 



remnant as well as a more complete consideration of the nu- 
clear composition of the ejecta. 
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